Method for frequency domain seismic data processing on a massively parallel computer

ABSTRACT

A frequency domain method of processing geophysical data on a computer having massively parallel processors. The method involves assigning data slice partitions to each processor, precomputing a velocity model corresponding to the geophysical data, and migrating the data on each slice within each processor using a one-pass, split wave equation finite difference technique for depth migration and either phase shift or recursive techniques for time migration. A sequence of transforms and transpositions between processors assigned partitions on each frequency slice transforms into the frequency-wavenumber domain and allows the migration calculations to be directly performed by each processor to be independent of each other processor. The transforms and transposes also allow for depth migration error correction and filtering in the frequency-wavenumber domain.

FIELD OF THE INVENTION

This invention relates to the field of geophysical prospecting.Specifically, the invention involves a frequency domain method ofprocessing seismic data using parallel processors.

BACKGROUND OF THE INVENTION

The search for subsurface hydrocarbon deposits typically involves asequence of data acquisition, analysis, and interpretation procedures.The data acquisition phase involves use of an energy source to generatesignals that propagate into the earth and reflect from varioussubsurface geologic structures. The reflected signals are recorded by amultitude of receivers on or near the surface of the earth, or in anoverlying body of water. The received signals, which are often referredto as seismic traces, consist of amplitudes of acoustic energy whichvary as a function of time, receiver position, and source position and,most importantly, vary as a function of the physical properties of thestructures from which the signals reflect. The data analyst uses thesetraces along with a geophysical model to develop an image of thesubsurface geologic structures.

The analysis phase involves procedures that vary depending on the natureof the geological structure being investigated, and on thecharacteristics of the dataset itself. In general, however, the purposeof a typical seismic data processing effort is to produce an image ofthe geologic structure from the recorded data. That image is developedusing theoretical and empirical models of the manner in which thesignals are transmitted into the earth, attenuated by the subsurfacestrata, and reflected from the geologic structures. The quality of thefinal product of the data processing sequence is heavily dependent onthe accuracy of these analysis procedures.

The final phase is the interpretation of the analytic results.Specifically, the interpreter's task is to assess the extent to whichsubsurface hydrocarbon deposits are present, thereby aiding suchdecisions as whether additional exploratory drilling is warranted orwhat an optimum hydrocarbon recovery scenario may be. In thatassessment, the interpretation of the image involves a variety ofdifferent efforts. For example, the interpreter often studies the imagedresults to obtain an understanding of the regional subsurface geology.This may involve marking main structural features, such as faults,synclines and anticlines. Thereafter, a preliminary contouring ofhorizons may be performed. A subsequent step of continuously trackinghorizons across the various vertical sections, with correlations of theinterpreted faults, may also occur. As is clearly understood in the art,the quality and accuracy of the results of the data analysis step of theseismic sequence have a significant impact on the accuracy andusefulness of the results of this interpretation phase.

In principle, the seismic image can be developed using athree-dimensional geophysical model of seismic wave propagation, therebyfacilitating accurate depth and azimuthal scaling of all reflections inthe data. Accurately specified reflections greatly simplify datainterpretation, since the interpretational focus can be on the nature ofthe geologic structure involved and not on the accuracy of the image.Unfortunately, three dimensional geophysical models frequently requireintolerably long computation times, and seismic analysts are forced tosimplify the data processing effort as much as possible to reduce theburdens of both analysis time and cost.

In addition to the 3-D computation challenge, the analyst faces aprocessing volume challenge. For example, a typical data acquisitionexercise may involve hundreds to hundreds of thousands of sourcelocations, with each source location having hundreds of receiverlocations. Because each source-receiver pair may make a valuablecontribution to the desired output image, the data handling load (i.e.,the input/output data transfer demand) can be a burden in itself,independent of the computation burden.

Seismic data analysts have historically used several differentapproaches to manage these burdens, directly or indirectly. Theseapproaches relate principally to either the manner in which the dataacquisition exercise is designed and carried out, or to the assumptionsmade during the data analysis effort. In both cases, the quality of theoutput of the data interpretation procedure may be directly affected.These approaches are most easily discussed in conjunction with FIG. 1,which depicts a perspective view of a region 20 of the earth for which ageophysical image is desired. On the surface 18 of the earth are shown anumber of shot lines 2 along which the seismic data are acquired. Asshown in FIG. 1A, shot lines 2 consist of a sequence of positions atwhich a seismic source 3 is placed and from which seismic signals 5 aretransmitted into the earth. Receivers 4 placed along each line receivethe signals from each source position after reflection from varioussubsurface reflectors 6.

A first method of managing the seismic data burdens discussed aboveinvolves careful definition of the region over which the data areacquired. Specifically, use of any available preliminary geologic andgeophysical information may facilitate the minimization of the surfacearea over which seismic data may need to be acquired. Such aminimization will directly reduce the amount of data that is ultimatelyacquired. Furthermore, similarly careful planning of the spacing betweenshot lines will optimize the analysis effort by reducing data volume.And finally, optimization of the number of sources and receivers thatare used, and of the spacing between adjacent source and receiverpositions, will also benefit the data analyst.

None of these efforts can be accomplished without a penalty. Forexample, relatively wide spacing between shot lines, or between sourcesand receivers, reduce the resolution of the computed seismic image, thusmaking interpretation more difficult. In addition, complex geologicfeatures may not be resolvable without relatively close spacing. Andfinally, certain data acquisition exercises, such as in relativelyunexplored areas, do not allow optimization of the surface area overwhich data is to be acquired. As a result, the data handling burdencannot be entirely eliminated through data acquisition planning.

Methods of minimizing the computational burden are often implementedduring data analysis. One commonly invoked technique involves use of atwo-dimensional geophysical model. For example, in FIG. 1A, the signalsfor each source are depicted as traveling in the plane directly beneaththe shotline on which the source lies. Thus, the signal is assumed topropagate independent of out-of-plane geologic structures. Thissimplifying assumption allows use of two dimensional geophysical modelsin the image generation process, and, as is well known, two dimensionalanalysis procedures can be much more computationally efficient thanthree dimensional analysis procedures.

Limitations to the 2-D analysis assumption exist. Geologic structuresare rarely, if ever, two dimensional; that assumption may therefore leadto inaccurately specified images. Because little is generally known ofthe geologic structure being investigated, the analyst usually does notknow the extent to which that image is in error. In addition, becauseeach plane is analyzed independently, the interpreter must tie theimages for each plane to each of the others by interpolation or othersimilar interpretative methods if a continuous image across the entirecubic region is desired. Finally, some complex structures, such asfaulted regions and salt features, cannot be accurately analyzed merelyby use of two dimensional methods.

Because of these and other limits that have long constrained seismicdata analysts, the petroleum industry has typically been an early userof newly developed high speed computer hardware. As each new generationof equipment has become available, analysis routines that implementfully three dimensional analysis capabilities have become more commonlyused. Nevertheless, it is not uncommon for significant computer times tobe involved in complex analyses, often involving weeks or months ofactual processing time.

The recent availability of massively parallel processors offers asignificant opportunity to seismic data analysts. Massively parallelprocessors (MPPs) have multiple central processing units (CPUs) whichcan perform simultaneous computations. By efficient use of these CPUs,the weeks or months previously required for complex analyses can bereduced to a few days, or perhaps a few hours. However, this significantadvantage can only be realized if efficient computational algorithms areencoded in the MPP software. Thus, the opportunity MPPs offer seismicdata analysts also creates a challenge for the development of suitablecomputational algorithms that take advantage of the multiple CPUs.

This challenge can be easily discussed by considering the manner inwhich computational algorithms have most commonly been written forexisting seismic analysis routines. Until recently, computers relied ona mode of operation referred to as sequential computing. Sequentialcomputing involves use of analytic routines that perform only a singleprocedure, or perhaps focus on a single subset of the data or image, atany given time. This is a direct result of a computer having only oneCPU. For that reason, the only optimization procedures that can beemployed on single CPU computers are those which increase the efficiencyof the processing as to the procedure or subset. Because allcalculations must ultimately be performed by that single CPU, however,the options for obtaining high performance are innately limited.

On the other hand, the multiple CPU capability of MPPs offers an obvioussimultaneous computation advantage. This advantage is that the totaltime required to solve a computational problem can be reduced bysubdividing the work to be done among the various CPUs, provided thatthe subdivision allows each CPU to perform useful work while the otherCPUs are also performing work. Unfortunately, the disadvantage ofmultiple CPU hardware is that the sequential processing methods thathave long been used in software development must be replaced by moreappropriate parallelized computing methods. Simply stated, MPPs requirethat processing methods be developed which make efficient use of themultiple CPU hardware. Ideally, these methods should organize thedistribution of work relatively evenly among the processors, and ensurethat all processors are performing necessary computations all of thetime, rather than awaiting intermediate results from other processors.

The challenge of defining parallelized processing methods, and ofoptimizing those parallelized methods once defined, is particularlyacute in the seismic data processing arena. Seismic data consistsgenerally of a large number of individual traces, each recorded somewhatindependently of the other traces. Logically enough, sequentialcomputing methods that require the analytic focus to be placed on asingle calculation at a time adapt well to analysis of these independenttraces. This is true even though computational bottlenecks may exist.For example, portions of the analytic sequence may require relativelymore computation time than other portions, must be completed beforeother calculations may proceed, or may rely on similar input data asother traces, for example traveltimes. Since no simultaneouscomputations occur in sequential processing, none of these bottleneckslead to a reduction in computational efficiency with a single CPU,except as to the total processing time that is required. Except as tothat total time requirement, the existence of such computationalbottlenecks does not otherwise pose problems for the analyst. To takefall advantage of MPP computing capabilities, however, where the goal isto perform simultaneous processing in all CPUs, methods for optimizingthe seismic analysis phase by eliminating such bottlenecks must bedeveloped.

This advantage of an MPP becomes clear by considering the limitationwhich calculation time places on image region size in single CPUcomputers. Increasing the size of the image, e.g., by expanding the sizeof cube 20 in FIG. 1, or increasing the amount of data to be processed,e.g., by adding additional sources 3 and receivers 4 to shotlines 2,increase the total computation. That direct impact on calculation timeplaces a heavy burden on seismic analysts to optimize image size,especially since even small image regions may require weeks ofcomputation time on even the highest speed sequential processingcomputers. In contrast, efficient processing on MPPs, which may have asmany as or more than 256 individual CPUs, should only involve minimallylengthened computation times, since each CPU would assume just afraction, for example 1/256, of the additional work required by thelarger region. This potential for scalability of the image region andthe work load required in image generation is a principal benefit ofMPPs, a benefit that can only be realized if parallelized seismicprocessing methods allowing such workload scalability are developed.

Basic considerations for determining efficient parallelized seismicprocessing methods become evident by reconsidering the above review ofthe seismic analysis process. As noted, the purpose of seismic analysisis to analyze measured seismic data using geophysical models to developimages of the subsurface. Therefore, each of three principal processingcomponents--data, model, and image--may be considered to be a candidatefor distributing computational work among the various processors in anMPP. One option for distributing work among the processors would be toassign different groups of the input seismic trace data to differentprocessors. For example, traces may be grouped by source locations, withdifferent processors being assigned different groups. Similarly, theoutput image could be subdivided and assigned to different processors.Finally, it may also be possible to subdivide the geophysical model usedto generate the output image into groupings that can be assigned to thevarious processors. (That model is generally considered to be embodiedin the arithmetic operations required by the mathematical model that isthe subject of the processing effort. For example, in seismic analysisthe mathematical model is often based on the wave equation). Forexample, the data may be transformed into the frequency domain, withindividual frequencies assigned to individual processors. It may also bepossible to develop combinations of these approaches. For example,groups of processors may be assigned collective responsibility forspecific frequencies in the model and all depths in the image, whilehaving individual responsibility for specific horizontal locations inthe image. The challenge to the seismic data analyst is to determinemethods of subdividing the seismic data, model, and image intocomponents that can be assigned to individual processors in the MPP,thus allowing calculations to be performed in each processorindependently of other processors. This subdivision of seismic dataanalysis into individual components is commonly referred to as seismicdecomposition.

One type of MPP has from thousands to tens of thousands of relativelyunsophisticated processing elements. In this kind of machine, theprocessing elements typically perform the same operation on multipledata streams, a Single Instruction, Multiple Data stream (SIMD) machine.An example is the CM2, a product of the Thinking Machines Corporation.These kinds of machines typically lack shared memory, i.e., eachprocessor has its own separate memory unit and the information in thememory cannot be directly accessed by other processors. The individualprocessors typically have limited computing capability and memory.Because of the large number of processing elements and a lack of sharedmemory, data transfer between the processing elements is a majorbottleneck in efficient utilization of the capability of the machines.Even with sophisticated interconnection techniques, such as in ahypercube arrangement, transfer of data between processors is a majorfactor in the running time of programs.

Other computers have much more powerful elements in arrays of tens orhundreds. The T3D, a product of Cray Research Corporation, is an exampleof this kind of machine. Besides having individual processing elementsthat are much more powerful than those in the CM2, the T3D has fewer ofthe elements and a physically distributed, logically shared memory. ThisMultiple Instruction Multiple Data stream (MIMD) machine has differentelements performing different operations on different parts of the dataat the same time. The reduced number of processing elements means thatdata does not have to be transferred to as many elements as in a SIMDmachine. Because of the increased sophistication and cost of theindividual elements and because of their fewer numbers, efficientutilization requires that the load on the processing elements bebalanced. An additional factor is that each processing element mustaccommodate a larger subset of the overall data volume; computationsthat involve sorting of the data could become more complicated.

U.S. Pat. No. 5,404,296 issued to Moorhead discloses and claims a methodfor migration on an MPP in which there are a large number of processingelements arranged in a preselected, regular pattern. The data areinitially shot-sorted and partitioned into blocks of shots so that theproduct of the number of receiver positions in the data set and thenumber of shots in a block equals the total number of processingelements available. Because 3-D seismic data volumes typically containhundreds of shot and receiver positions, the disclosure and claims arelimited to SIMD machines.

One approach that reduces the amount of the 3-D processing is to image apartially processed data set. This is particularly useful when thepartial processing operation is computationally simpler than a completeimaging process. For example, prior to migration, the data could bestacked using a conventional Normal Moveout (NMO) velocity analysis. Aswill be familiar to those knowledgeable in the art, algorithms thatperform migration of zero-offset data are computationally less intensivethan those that migrate unstacked data. The stacking operation generallyimproves the signal to noise ratio of the data compared to that of theindividual traces. To the extent that the stacked data accuratelyrepresents a hypothetical zero-offset trace, the computational burden isreduced.

In U.S. Pat. No. 5,349,527, Pieprzak and Highnam disclose a method formigration of such a data set on a SIMD MPP. The method starts with a 3-Ddata volume consisting of a time series (vertical axis) of reflectionseismic data on a regular rectangular grid (horizontal axes) on thesurface of the earth [(t, x, y) data domain]. A Fourier transformationis performed on the reflection seismic data to give a 3-D volumeconsisting of amplitudes as a function of frequency (vertical axis) onthe rectangular surface grid [(ω, x, y) data domain, where ω is thetemporal frequency]. The data are partitioned into frequency sub-bands,and within each frequency sub-band, at selected frequencies, data arepartitioned into subsets of the rectangular surface grid and assigned toindividual processors. This is necessitated by the limited memorycapability of the individual processors. The limitations of theindividual processors also necessitate a hypercube arrangement of theprocessors in order to reduce the time involved in data transfer betweenthe processors.

In Pieprzak and Highnam's invention, the subsequent migration isaccomplished by an iterative, two-step process. The first step is thedownward continuation of a single frequency component over one depthinterval. The contributions of all the frequencies are then summed toproduce a migrated image at the bottom of the depth interval. The twosteps are repeated for successive depth intervals until the maximumdepth is reached.

The downward continuation is performed using the McClellan transformmethod given in Hale, D., 3-D Depth Migration via McClellanTransformations, 56 Geophysics 1778 (1991). The actual implementation isdone by means of a convolution filter. As will be known to thosefamiliar with the art, the Hale method has difficulty handling steepdips and the computational burden becomes very heavy at dips greaterthan 60°. To reduce the computational burden, Pieprzak and Highnam teachthe use of a recursive Chebyshev filter for the convolution step. Thecoefficients of the filter are precomputed from the velocity model ofthe subsurface and stored redundantly in the local memory of theprocessing elements. Such an approach is appropriate in SIMD computerswhere the processing elements need to keep track of only small subgridof (x, y) points and do not have the capability of computing the filtercoefficients "on the fly."

Those knowledgeable in the art recognize that for large problems,frequency domain methods can offer significant computation cost savings.Furthermore, the prior art two-pass migration methods arecomputationally fast because they perform the migration in the x- andy-direction separately. Pieprzak and Highnam use a one-passfrequency-space domain migration process that is computationally moreaccurate than the two-pass migration process used in prior art. However,the increased accuracy of the one-pass method comes at the cost ofvastly increased computational load. One option for overcoming thatincreased computational load is to employ frequency-wavenumbermigration, which can be significantly faster than frequency-spacemigration. However, a frequency-wavenumber scheme based on Pieprzak andHighnam's disclosure would be impractical for the reason that thetransformation from the spatial domain to the wavenumber domain becomesexpensive when each processor has access to only a small portion of thex-y grid.

U.S. Statutory Invention Registration H 482 issued to Berryhill,Gonzalez and Kim in 1988 discloses a method for recursive time migrationin the frequency-wavenumber domain. That method assumes that thesubsurface of the earth can be modeled by layers in which the velocityis constant. However, the method is inefficient to implement on an MPPbecause it requires successive Fourier transformations in differentdirections on a data volume entirely within processor memory, and thegeophysical data sets of interest to industry cannot be stored entirelywithin processor memory.

It would therefore be desirable to have a method that is able to performa one-pass migration in the frequency-wavenumber domain on an MPP in acost-effective manner. The present invention satisfies this need.

SUMMARY OF THE INVENTION

The present invention is a method of processing seismic data on parallelprocessors, preferably on a massively parallel processor (MPP). The endproduct of the invention is a migrated seismic section in the form of animage of the subsurface of the earth. The input to the process may be ofeither stacked seismic data or of Dip Moveout (DMO) corrected commonoffset seismic data. The process migrates seismic data in thefrequency-wavenumber domain, thereby significantly increasing thecomputational speed over frequency-space methods.

A velocity model is used by the individual processors of the MPP toobtain parameters for the migration operation. The individual processorsperform migration calculations independently of the other processors.The invention partitions the data on the individual processors to makethe subsequent processing efficient. The process uses a series of matrixtranspositions to efficiently implement the Fourier transformationsrequired for the migration in an efficient manner on an MPP.

In one embodiment, the data are migrated to produce an image, in depth,of a 3-D data volume. This embodiment uses a computational scheme inwhich the velocity of propagation of seismic waves in the subsurface canvary laterally and vertically. The data are downward continued in thefrequency-space domain from one depth to the next using the velocitymodel; summation of the various frequency components produces a migrateddepth image. The downward continuation is done in two steps: First, aparaxial approximation to the wave equation is solved using an implicitfinite difference method with a splitting of the equation in the inlinex and crossline y directions. This is analogous to prior art that used atwo-pass migration process. The splitting reduces the computationalburden by several orders of magnitude but introduces an error in thesolution. Second, an error correction is performed that compensates forthe error introduced by the process of splitting the wave equation intox and y components. The error correction is computationally fast becauseit is performed in the frequency-wavenumber domain, and does not have tobe carried out at every depth interval. The result is an accuratesolution of the wave equation valid for steep dips.

The solution of the wave equation for the downward continuation steprequires the solution of a large tridiagonal system of equations. Thepresent invention includes a novel "burn at both ends" method forsolution of such a system of equations on an MPP. By efficient use ofthe processors, the method performs the downward continuation inapproximately half the computation time that would be needed in priorart methods of solving such a system of equations.

In a second embodiment, the end product of the process is a migratedtime image. The seismic data (stacked traces or DMO corrected traces)for a 3-D volume, after pre-processing, are downward continued on an MPPcomputer. The method starts at the surface and proceeds time-step bytime-step using a phase shift method. Summation of the contributions ofthe individual frequencies produces a migrated image. The secondembodiment is implemented in the frequency-wavenumber domain, therebyproviding a significant increase in computational speed over prior art.

In a third embodiment, the end product of the process is also a migratedtime section. The velocity model is a piecewise continuous function ofdepth. Starting at the surface with data in a (t, x, y) cube, theso-called Stolt migration is applied over each interval. This is donerecursively, alternating the Stolt migration with a downwardcontinuation of the wavefield in each layer. This embodiment is alsoimplemented in the frequency-wavenumber domain. For subsurface modelsthat have relatively thick layers of constant velocity, this embodimentoffers a significant improvement over the second embodiment.

BRIEF DESCRIPTION OF THE DRAWINGS

The present invention and its advantages will be more easily understoodby reference to the following detailed description and the attacheddrawings and tables in which:

FIGS. 1 and 1A schematically illustrate the data acquisitionconfiguration typically used in 2-D seismic exploration.

FIG. 2 shows the preprocessing of the data required for a depthmigration embodiment.

FIGS. 3A and 3B show the overall flow of data in the invention andconceptual steps involved in a first depth migration embodiment.

FIG. 4 illustrates the steps involved in the downward continuation andimaging calculations for the embodiment depicted in FIGS. 3A and 3B.

FIG. 5A, 5B and 5C pictorially depict the broadcasting of velocityslices between processor groups, the partitioning of a frequency sliceamong the PEs in a to different processor group and the downwardcontinuation of a single frequency slice from the surface down one depthinterval.

FIG. 6A and 6B show the procedure for retrieving previously computedpartial images and the method for combining partial images in a depthmigration logarithmic fan-in.

FIGS. 7A, 7B, 7C, 7D and 7E show the overall flow of data in timemigration embodiments of the present invention.

FIG. 8 depicts the calculation time improvements for a depth migrationembodiment of this invention as compared to a prior art technique.

DETAILED DESCRIPTION OF THE INVENTION

The present invention is a method of processing seismic data on aparallel processor. The method may be employed on any type of parallelprocessing computer that has more than one processor available toperform analysis and/or input/output tasks. For example, a Cray Y-MP,which has up to eight processors may be employed. Networked personalcomputers or workstations may also be employed. In a preferredembodiment, the method will be employed on a massively parallelprocessor (MPP) having 64 or more processing elements. One example ofsuch an MPP is the T3D made by the Cray Research Corporation. Forconvenience, and not to be construed as limiting, the abbreviation MPPwill be used herein to refer to any parallel processing computersuitable for the present method. However, use of MPP in that manner isin no way to be construed as limiting application of the present methodto solely those computers deemed by the commercial marketplace to bemassively parallel.

While the implementation is discussed in terms of a MIMD machine, theinvention is not restricted to them. Apart from velocity data, thepresent method involves use of a relatively small number of inputprocessing parameters, and therefore can also be usefully implemented onSIMD MPPs. In both MIMD and SIMD implementations, the input parametersare read into one processor and broadcast to the other processors.

As is well known, seismic data processing methods generally start withdata from a portion of the earth's surface. FIG. 1 depicts a perspectiveview of a region 20 of the earth for which a geophysical image isdesired. On surface 18 of the earth are shown a number of shot lines 2along which the seismic data are acquired. As shown in FIG. 1A, shotlines 2 consist of a sequence of positions at which a seismic source 3is placed and from which seismic signals 5 are transmitted into theearth. Receivers 4 placed along each line receive the signals from eachsource position after reflection from various subsurface reflectors 6.As will be known to those familiar with the art, it is now common forthe acquisition to be carried out with multiple rows of sources andreceivers. The present disclosure is meant to encompass these situationsas well and the preliminary processing steps involved in the "binning"of such a data set are part of the art.

These reflection signals will typically be preliminarily processed usingone of two known techniques. First, the preliminary processing may focuson producing a stacked data volume. The process of stacking consists oftaking data from single line 2. By a process of velocity analysis andnormal moveout (NMO) correction, an approximate zero-offset seismicsection is obtained. In this 2-D section, the horizontal (x) coordinateis the position along the line 2, the vertical coordinate is time, andthe data represent reflection amplitudes. By similar processing of anumber of parallel lines, a 3-D data volume is obtained.

A second preliminary processing alternative may be used in which thedata along the seismic lines may be dip moveout (DMO) corrected toproduce a series of common offset panels with the same shot-receiverdistance. The end result of both stacking and DMO correction, as is wellknown to those skilled in the art, is a 3-D data volume of reflectionamplitudes in (x, y, t) space.

As would be familiar to those knowledgeable in the art, the processingdirection need not be the same as the direction of data acquisition.Furthermore, it is possible to acquire seismic data from a single sourcewith a multiplicity of lines of receivers. Such an acquisition geometrymakes it possible to obtain a 3-D image of the earth and to perform thestacking and DMO correction in arbitrary directions. The presentdisclosure is intended to include such datasets.

For convenience, the detailed description of several embodiments of thepresent invention will be organized as follows: First, an introductorysection discusses preprocessing considerations. Second, an overview ofthe major processing steps for the embodiment is presented. Second, adetailed description of the data manipulation procedures necessary forefficient implementation of the embodiment on an MPP is presented.Finally, the imaging scheme used in the embodiment is discussed. Thisorganization is not be construed as limiting, however, but is merelyused to simplify discussion of the various aspects of the presentinvention.

Depth Migration

Depth migration embodiments of the present invention allow for bothefficient implementation of the depth migration imaging scheme in thefrequency-space domain and for the optional implementation of an errorcompensation scheme to correct errors introduced by the migration. Theerror compensation is performed in the frequency-wavenumber domain,thereby allowing accurate migration results to be obtained in reducedMPP computation times. The embodiments include a novel solution schemefor the innermost portion of the migration calculation which reducescalculation times as compared to the prior art.

Preprocessing Considerations

As noted above, the preliminary processing of the seismic data depictedin FIG. 1 will generally involve either stacking or DMO analysis. Theresult of that preliminary processing is a 3-D data volume consisting oftime series of reflection amplitudes recorded on a surface grid ofinline x and crossline y directions, thereby defining the (t, x, y) datavolume 30 depicted in FIG. 2. Data volume 30 consists of time series(vertical axis) recorded at various points on a horizontal surface grid,32. Depth migration embodiments of the invention will generally continuewith an additional data pre-processing step, performed on a computersuch as the Cray Y-MP, to transform the data into a 3-D volume in (x, y,ω). This step involves a Fourier transformation into a volume 34 inwhich the vertical axis is frequency and the horizontal axes areunchanged. This step will be understood to those skilled in the art.

Next, volume 34 is sliced perpendicular to the ω axis to give frequencyslices, 36, to be used in the migration. The number of slices which aregenerated and the frequencies for each slice will be a function of theanalysis to be performed, the nature of the dataset being analyzed, andother factors which are well known. The slices may optionally becompressed to minimize memory and data input/output communicationburdens. The data compression might consist, without limitation, of anynumber of well known procedures, such as representing the data withfewer bits.

Depth migration preprocessing may also slice a precomputed depth modelfor the migration. The depth model will consist of velocities as afunction of depth on the same (x, y) grid as the data volume, and thedepths and number of velocity slices will be also determined by theanalyst. The velocity slices, 38 in FIG. 2, may also be compressed. Thefrequency slices and velocity slices are written to disk or other highspeed device for subsequent use by the MPP. Hereinafter, the term "disk"is intended to include other high speed devices on which data can bewritten and from which it can be retrieved. Each of the above stepswould be familiar to those knowledgeable in the art.

An initial allocation of the frequency slices to the processors in theMPP is performed based upon the number of x and y points in the data,the number of frequencies to be used in the migration and the memoryavailable on each processor in the MPP. For convenience, this allocationis performed as part of the startup routine assigned to one PE withinthe MPP which acts as a control processor for the MPP in determining theallocation. The allocation will typically involve determining theminimum number of processors required per frequency slice, with eachprocessor having an equal number of x and y data points at whichcalculations will be performed. However, the present method is notlimited to use solely of that minimum number or to maintaining thatequal number per processor, which are selected herein for convenience.Upon completion of this allocation, the control processor returns tofunctioning as an analysis processor, in the same manner as all otherprocessors in the MPP. The purpose of the allocation is to determine thenumber of processors which will be required to perform computations foreach frequency slice. The result of the allocation is that, in thecomputational procedures and data manipulations discussed further below,each processor is assigned one or more frequency slices, or in thealternative, a group of processors will be assigned a single frequencyslice. If a single processor can be assigned more than one frequencyslice, the frequencies are chosen so as to balance the computationalload on the processors. Ordinarily, this will mean that a low frequencyand a high frequency slice will be assigned in pairs to the processor.As would be familiar to those knowledgeable in the art, migration of lowfrequencies is less of a computational burden than migration of highfrequencies; therefore such a pairing helps balance the computationalload among the processors.

Processing Overview

Data Flow--Generally

Referring to FIG. 3A, at "START" the data are read from disk to the PEs.If the number of frequency slices per processor, 42, is less than one,the processing follows the steps in FIG. 3B. If that number is one ormore, then the frequency slices are collected into sets, 44, areassigned to processing elements (PEs), and read into the memory of eachPE, 46. If the PEs are processing an initial set of frequencies in afirst pass analysis, 48, the process of downward continuation andimaging (to be described below in conjunction with FIG. 4) is applied,50. The result of this downward continuation and imaging is to create apartial image plane for each depth.

If additional frequencies remain to be processed, as is usually thecase, additional frequency slices are assigned and read into each PEfrom disk, 46. For the processing of the second and subsequent set offrequencies, 48, the process of downward continuation and imaging, 54,updates the previously created image plane for each depth. This sequencerepeats as long as additional frequencies are to be processed, 52. Themigration is complete when all frequencies have been processed, 55.

If more than one processor is required to hold the data for a singlefrequency slice, the steps in FIG. 3B are followed. FIG. 3B is analogousto FIG. 3A except that groups of processors are required to process eachfrequency. Each frequency slice is assigned to a group of PEs and readinto each PE in the group, 56. For the first set of frequency slices,58, the process of downward continuation and imaging (to be describedlater in conjunction with FIG. 4) is applied, 60. As above, the resultof this downward continuation and imaging is to create an image planefor each depth. If additional frequencies are to be processed, as isusually the case, these additional frequencies are assigned and readinto each group of PEs, 56. Also as above, for the processing of thesecond and subsequent set of frequency slices, 58, the process ofdownward continuation and imaging, 64, updates the previously createdimage plane for each depth. If all the frequencies have been processed,62, the process is terminated, 65.

One of the noteworthy features of this invention is that depth migrationembodiments have the ability to store intermediate results after a setof frequencies have been processed, 55 in FIG. 3A and 65 in FIG. 3B. Themigration can therefore be continued or restarted at a later time,allowing correction or refinement of input parameters and analysis ofintermediate results. The desirability of such a restart feature will befamiliar to those knowledgeable in the art, although the cost effectiveimplementation of such a feature has not previously been attained.

Data Flow--Downward Continuation and Imaging

The downward continuation and imaging step, (50 and 54 in FIG. 3A, 60and 64 in FIG. 3B), is described in more detail in FIG. 4. This processis performed by each processor or group of processors, i.e. 46 in FIG.3A and 56 in FIG. 3B. The process starts, 70, at an initial depth in thesubsurface, usually the surface of the earth. Next, the velocity datafor the underlying depth interval is read from disk into one of theprocessor groups, 72. (It is to be understood from here on that whenreference is made to processor groups, a group may consist of just oneprocessor when one processor holds one or more frequencies.)

Because all processor groups in the MPP will be performing the downwardcontinuation, efficient utilization of the MPP requires that thisvelocity information be efficiently made available to all the processorgroups. In the preferred embodiment, the broadcast of velocities, 74, isperformed by a process of logarithmic fanout. This step of logarithmicfanout is illustrated in FIG. 5A. Logarithmic fanout enables thevelocity information needed for the downward continuation to be rapidlydisseminated to all the processors. A compressed velocity slice, 38,from disk is loaded into a group of processors, labeled #1, anduncompressed. Subsequently, the uncompressed velocity data from #1 aretransferred to a second processor group, #2. At the next stage, #1 sendsthe velocity slice to #3 while simultaneously, #2 sends it to #4. At thenext stage, #1 sends velocity information to #5, #2 sends it to #6, #3to #7 and #4 to #8. In this fashion, the number of processors that havethe velocity information grows geometrically. This is known as alogarithmic fanout.

The labeling of the processor groups in FIG. 5A is for illustrativepurposes only. In reality, the transfer of velocity data may be betweenprocessors in the same physical neighborhood within the MPP, or whichmay be otherwise closely connected by the MPP's internal communicationpathways, or based on another communication relationship applicable tothe MPP system being used. It is also to be understood that otherschemes for dissemination of the velocity data could be used and theinvention is not restricted to a logarithmic fanout. These would befamiliar to those knowledgeable in the art.

In some MPPs, the processing elements have a buffered read and writecapability. This embodiment of the present invention takes advantage ofthis capability by reading the velocity data from disk asynchronously.While the velocity data for one depth are being broadcast to all theprocessors, the velocity data for the next depth interval are being readfrom disk into the first processor. This procedure would be familiar tothose knowledgeable in the art.

At the completion of step 74 of FIG. 4, each PE will have frequencyslice data in its memory, via the assignment and reading procedures ofstep 46 in FIG. 3A and step 56 in FIG. 3B. Each PE will also have thevelocity data from step 74 of FIG. 4. If the number of frequencies perPE is greater than or equal to one, then each PE will have the frequencyand velocity slice data for all x and y positions in the data volume inits memory. However, if more than one processor is needed for afrequency slice, each PE will be assigned, and will retain in itsmemory, frequency slice and velocity data for the full range of the datavolume along the x axis, but for only a limited range of y values. Forexample, a single frequency slice, 36 from FIG. 2, is shown in FIG. 5B.The partitioning of slice 36 is shown in 36', where each partitionextends the entire length of the x-axis, but less than the entire lengthof the y-axis.

The downward continuation step, 76 in FIG. 4, consists of the PEs foreach frequency slice 36 performing the downward continuation calculationfor its frequency ω_(n). The downward continuation process, 76, takesthese data at the initial depth, combines it with the velocityinformation in the underlying depth interval and produces frequency datafor the bottom of the depth interval. FIG. 5C depicts the downwardcontinuation from the surface down one depth interval. Frequency slice,36', corresponding to frequency ω_(n), is downward continued from depthz=0 (92) to depth z=Δz (94). Each PE will perform the downwardcontinuation for an entire frequency slice, or for the entire partitionit has been assigned if more than one PE is required for each slice.

The downward continuation procedure in step 76 may involve any number ofknown procedures, as is discussed further below. In one embodiment, anapproximation to the wave equation is used which reduces thecomputational complexity of the downward continuation process. Theapproximation relies on a paraxial form of the wave equation and asplitting in the inline x and crossline y directions. This is followedby an error correction scheme, 78 in FIG. 4. The error correction schemecompensates for errors introduced by the paraxial approximation, by griddispersion and by the splitting of the wave equation into x and ydirections. The error correction will not generally need to be done atevery depth for which a downward continuation is performed, depending onthe velocity model's characteristics. The downward continuation anderror correction processes are discussed further below.

As will be understood by those skilled in the art, the summation, 80, ofthe contributions from the various frequency slices produces a depthimage at the bottom of this depth interval. Where prior images existfrom other frequencies, (as a result of prior analysis in 50 and 54 inFIG. 3A, or in 60 and 64 in FIG. 3B), these prior images are combinedwith the image produced at this step, 80. Steps 76, 78, and 80 aresometimes referred to by those skilled in the art as the migrationprocess.

The process of summation of the frequency components and the priorimages is illustrated in FIGS. 6A and 6B. FIG. 6A shows a series ofpartial images, 100, in the data volume that have been computed inearlier steps for depths Z₁ -Z_(n). The partial image corresponding tothe depth at which the most recent calculations have been completed,shown as 100' in FIG. 6A, is retrieved from disk and transferred to aPE, labeled in FIG. 6A as #1. As noted above, PE #1 may also be a groupof PEs. FIG. 6B shows PE #1 and other PEs, labeled #2 through #8, whichhave most recently performed downward continuation for a set offrequencies ω₁ through ω₈. The previously existing partial image at thedepth of interest, 100', is added to the computed image 110 in PE #1 toproduce an update partial image, 110". Similarly, the downwardcontinuation images for frequencies ω₂ through ω₈ on processor groups #2through #8 (all of which are shown as 110 in FIG. 6B for convenience)are combined with partial image 110" in a scheme that is the reverse ofthe FIG. 5A logarithmic fanout scheme used in broadcasting velocities,with the result that a new image estimate 100'" is obtained. Imageestimate 100'" in FIG. 6B is an updated estimate of image 100' from FIG.6A at the depth corresponding to the depth of image 100'. Where morethan one frequency slice fits on a single processor, the images fromthese frequencies are combined within the processor prior to beingcombined with the images on other processors.

The updated image is written to disk or other hardware device, 82 inFIG. 4. It should be understood here that the number and labeling of theprocessors and frequencies are for illustrative purposes. The presentmethod may be implemented with image updates occurring among any numberof processors, and the reverse fanout scheme may be adapted accordingly.Other schemes may also be used to generate image updates. Finally, acheck is made to see if the maximum desired depth has been imaged, 84.If not, the process is repeated for successive depth intervals, 86, 72.

As part of the load balancing on the processors, the same processorgroup does not perform the tasks of: (i) reading in and broadcasting avelocity slice as shown in FIG. 5A, (ii) reading in a previouslycomputed image slice as shown in FIG. 6A, and, (iii) writing out animage slice to disk as shown in FIG. 6B. Instead, these tasks areperformed by different processor groups, i.e., the processor groupslabeled as #1 in FIGS. 5A, 6A and 6B may be physically different. Onereason for these different assignments is to increase the efficiency inthe use of processor memory. Specifically, in the present method, it ismore efficient if each processor has only one portion of its memoryallocated to be a data input/output buffer. Because that buffer can onlybe accessed for one task at a time, it is more efficient if eachprocessor is assigned either an input or an output role at a time, andtherefore the load balancing scheme of the present method involvesensuring that each processor is assigned only one I/O task at a time.

Data Manipulation

A well-understood challenge to the seismic data analyst is to developdata manipulation and storage schemes which efficiently andcost-effectively implement the processing procedures of interest on thecomputer system available to the analyst. The present invention includesa unique series of data manipulations which take advantage of thecapabilities of the MPP, and which reduce the cost and time required forcompleting an analysis as compared to prior art techniques. These datamanipulations are most easily discussed in conjunction with Table I,which shows the sequence of steps in which transform and transpositionoperations are performed by the present invention on the arrangement ofthe data volume in PE memory during the depth migration. Table I alsoshows the sequence in which the data are stored in processor memory atthe end of each operation. These operations correspond generally to 76,78, and 80 in FIG. 4, as discussed above.

The frequency slice data are stored in PE memory in (x, y, ω) matrixformat, which in Table I for convenience are referred to generally asthe 1, 2, and 3 axes. Because each frequency slice will involve a fixedfrequency on, this embodiment of the present invention involves asequence of transform and transposition operations relating solely tothe 1 and 2 axes.

The initial sorting of data, shown in step 1 of Table I, results fromthe Fourier transform of the data volume and the specification offrequency slices discussed above in conjunction with FIG. 2. In step 2,the downward continuation procedure is carried out, also as discussedabove. This procedure occurs for all x and y grid points stored in eachPE. If the frequency slice is partitioned over more than one PE, theprocedure is repeated for all x and y grid points in the partition inthe PE. Step 2 in Table I corresponds to 76 in FIG. 4.

The next analysis procedures to be carried out, the error correction andwraparound filter procedures, correspond to 78 in FIG. 8. Theseprocedures, steps 3 through 9 of Table I are optional to the extent thatthey are not performed at every depth, as discussed below. Theseprocedures are most efficiently carried out in the wavenumber domain,and therefore are preceded by a series of matrix transpositions andFourier transformations, as shown in steps 3, 4, and 5 of Table I. Thetranspositions are required because more efficient use of each processorin the MPP results if all Fourier transformations are performed alongsolely the axis--1, 2, or 3 in Table 1--in which the PE has the entirerange of data. In the present method, all such transformations occuralong the 1-axis, as follows. As discussed above in conjunction withFIG. 5B, each processor in a group has the entire range of x data in afrequency slice, but only has a portion of they data (i.e. all dataalong the 1-axis of Table I, but only a portion of the data along the2-axis). Each processor can therefore perform the step 3 Fouriertransform of the x-data independently of the other processors. Theoutput of that step is (k_(x),y, ω) format data. Next, in step 4, amatrix transposition of the 1 and 2 axes, which relies on MPP librarymatrix transposition routines, is performed. This transposition occursamong all processors in the group with the result that each processorhas an entire range of y data in its memory, but for only a portion ofthe k_(x) axis, and thereby permits the application of the FFT to y datain each processor in step 5. The result from step 5 is that both the xdata and the y data are in the wavenumber domain. Step 6 is theapplication of the error correction and an optional wraparound filter,which does not involve a transformation of the (k_(x) k_(y), ω) data.

Use of the splitting introduces errors into the wavefield extrapolation,as is described in Brown, D. L., Applications of operator separation inexploration seismology, 48 Geophysics 288 (1983), and would be familiarto those knowledgeable in the art. As shown by Li, Zhiming, Compensatingfinite difference errors in 3-D modeling and migration, 56 Geophysics1650 (1991), the wave equation splitting error correction can be handledto first order by an extra phase shift term in the wavefieldextrapolation. This error correction step (FIG. 4, step 78) requiresthat the data be transformed from the (x, y, ω) domain to thefrequency-wavenumber domain, therefore requiring steps 3, 4, and 5 ofTable I. This error correction step usually needs to be applied onlyonce every few depth extrapolation steps, as will be understood to thoseskilled in the art. Where there are large vertical velocity gradients,the error correction may have to be applied more frequently. It is to beunderstood that higher order corrections, though not used in thisembodiment, could be implemented where there are large velocitygradients in the horizontal direction.

After the error correction is applied, an optional wraparound filter mayalso be applied. This wraparound filter would be familiar to thoseversed in the art.

After the error correction and wraparound filter, inverse procedures areperformed to re-obtain the (x, y, ω) domain. The inverse FFT in step 7puts the data in (y, k_(x), ω) space. The inverse FFT algorithm is alsoapplied along the 1-axis. Thereafter, the step 8 transpositionfacilitates the step 9 inverse FFT which brings the data into the (x, y,ω) space. The summation over ω in step 10 produces the depth-migratedimage in (x, y, z) space for a single frequency, as discussed above inconjunction with step 80 of FIG. 4.

An important advantage of the present invention is that the above datamanipulation scheme allows each PE to compute a direct "on the fly"solution of the wave equation in the downward continuation step, asopposed to prior art methods in which the PEs computed an indirect,look-up table solution by reliance on a precomputed filter coefficienttable. The prior art methods were forced to perform that indirectsolution due to the high computational burden which was previouslyrequired to directly solve the wave equation. The present scheme botheliminates that high burden and avoids those precomputed tables, thusproviding more accurate processing results than have previously beenattained in seismic data processing methods. That accuracy isindependent of the degree of dip, thereby avoiding a principal concernof the prior art methods. In addition, by avoiding the pre-computationand storing of coefficient tables, the present scheme also allows ahigher percentage of processor memory to be dedicated to computationaltasks, i.e. by allowing larger and/or more partitions to be stored ineach processor, thereby facilitating increased computational throughputcompared to the prior art methods

                  TABLE I                                                         ______________________________________                                        Parallelization of Depth Migration Data Manipulation Tasks                    Processors are in the 1-2 Plane.                                                                        Axis                                                Step    Task             1        2   3                                       ______________________________________                                        1       Start            x        y   ω                                 2       Downward Continue in x & y                                                                     x        y   ω                                 3       FFT              k.sub.x  y   ω                                 4       Transpose        y        k.sub.x                                                                           ω                                 5       FFT              k.sub.y  k.sub.x                                                                           ω                                 6       Apply Error Correction and                                                                     k.sub.y  k.sub.x                                                                           ω                                         Wraparound Filter                                                     7       Inverse FFT      y        k.sub.x                                                                           ω                                 8       Transpose        k.sub.x  y   ω                                 9       Inverse FFT      x        y   ω                                 10      Sum over ω x        y   z                                       ______________________________________                                    

Image Calculations

The data are downward continued (FIG. 4, step 76) by using an implicitfinite difference scheme to solve the paraxial wave equation, forexample as discussed in, Claerbout, J. F., Imaging the earth's interior,Scientific Publications, Inc. (1985). As would be known to thosefamiliar with the art, a direct implicit solution of the finitedifference wave equation requires the solution of a set of (n_(x) •n_(x)×n_(x) •n_(y)) equations, where n_(x) and n_(y) are the number of pointsin the x and y directions of the data volume. In contrast, an explicitfinite difference scheme requires fewer computations but is numericallyunstable. The present invention obtains an approximate solution to thewave equation using an implicit scheme and by splitting the waveequation in the x and y directions. This scheme requires the successivesolution of a much smaller set of matrix equations of order n^(x) ×n_(x)and n_(y) ×n_(y). The difference between implicit and explicit finitedifference solution schemes will be understood to those skilled in theart.

Both sets of matrix equations involved in the downward continuationprocess require the solution of equations of the form

    A.sub.x =b                                                 (1),

where the matrix A is tridiagonal. The vector x here refers to genericunknowns in the equation, not solely to the spatial coordinate x in thedata volume. As is well known, the conventional solution of thisequation proceeds by starting with the first equation and successivelyperforming a series of eliminations to get to the last equation, andthen back-substituting the solution obtained at that point in a reversesequence to solve all the equations. This conventional solution workssatisfactorily on an MPP where all the data for a single frequency sliceare on a single processor because each of the processors can workindependently of the other processors. This is also the case for thesolution of the x-equations because, as described above, all the dataare available in order within the memory of each processor. However,where more than one processor is needed for a single frequency slice,i.e. an individual processor is responsible for only a portion of they-axis, as in 36' of FIG. 5B, direct application of the conventionalmethod on an MPP is inefficient. Specifically, a computationalbottleneck is created because only one PE can be solving the y-equationsfor each ω-slice, although each PE can solve all the x-equations. As aresult, whenever a frequency slice is assigned to a group of more thanone PE, the PEs in the group cannot be efficiently utilized.

The present invention is an improvement over prior art in that they-equations may be solved so that two PEs are nearly always active oneach ω-slice. Specifically, while the PE assigned to the smallesty-values begins eliminating the wavefield at those values, the PEassigned to the largest y-values also begins eliminating the wavefield,starting at the largest y-value and proceeding toward smaller y-values.As will be understood to those skilled in the art, this process involvesthe sequential elimination of unknown wavefield parameters. As each PEreaches the limit of its range of assigned y-values, it passesinformation describing the modified equations to the PE responsible forthe adjacent range of y-values. This process continues until thewavefield can be completely determined at an interior value of y. Aswill also be understood to those skilled in the art, at this interiorvalue the set of equations characterizing the wavefield will have beenreduced from an indeterminate set to a determinate set, thus allowingthe wavefield at this location to be determined. Thereafter, thewavefield at this value of y is used in a reverse elimination todetermine the wavefield at both larger and smaller values of y, againpermitting two PEs to be active at the same time. A special case ariseswhen an odd number of PEs span the range of y-values. In this instance,the PE responsible for the central y-partition is active by itself for aportion of the calculation. This "burn at both ends" dual solutiontechnique, as just described, avoids the transposition of the ω-slicewhich would otherwise be necessary to place the entire y-range withineach PE as well as the corresponding inverse transpose. The "burn atboth ends" technique could also be used in fine-grained MPPs in both thex and the y directions to reduce computation time. The computer code forthis method is given in the Appendix.

In summary, the end product of the depth migration process is anaccurate migrated image at considerably lower cost than can be obtainedby a simple frequency-space migration. This is a result of using anapproximate solution to the wave equation, utilizing the x-y splitting,using the "burn at both ends" method for solving a set of tridiagonalequations, and applying an error correction in the (k_(y), k_(x), ω)domain, each of which are efficiently implemented on an MPP.

Time Migration

Second and third embodiments of the present method perform timemigration on an MPP. The second embodiment performs a time migration inthe frequency-wavenumber domain, an improvement over prior arttime-migrations which are limited to the frequency-space domain. A thirdembodiment offers a significant improvement in computational speed overthe second embodiment when the velocity in the earth is constant overrelatively thick intervals. The difference between these two embodimentslies in the imaging calculations which are performed.

Preprocessing Considerations

In a manner similar to the depth migration embodiment discussed above,an initial determination is made on the number of processors per group.There are differences, however, in the manner in which the data slicesare determined, and in the data flow. Moreover, at different stages ofthe processing, the number of processors in a group may be different.These embodiments are discussed further below, using a three partbreakdown analogous to that used above. In the following discussion,however, the data manipulation procedures which are involved in timemigration will be discussed first, followed by the discussion of thedata flow and imaging calculation procedures. This sequence simplifiesthe discussion of the more complex series of transpositions andtransformations that are involved in the time migration embodiments, ascompared to the depth migration embodiment discussed above.

Data Manipulation An important part of the time migration embodiments isthe manipulation of the data and the manner in which it resides in thememory of the MPP. Both embodiments start with data in a (t, x, y) datavolume, such as is depicted in FIG. 2. The MPP time migration sequenceinvolves a series of transpositions and transformations analogous tothose performed in the depth migration embodiment. However, thetransposition and transformation sequence expands upon that which wasrequired for the depth migration embodiment because Fouriertransformations have to be made in both temporal and two spatialdomains, as compared to the depth migration transformations which wereonly in the two spatial domains. As discussed above, the purpose of thetranspositions is to ensure that all transformations are performed alongthe 1-axis, as that axis was defined above in conjunction with Table Iand as used below in conjunction with Table II. The desire to performall transformations along the same axis derives from the preference thatonly one PE be involved in the transformation, as was also noted abovein the depth migration embodiment. Table II shows the effects of thetransform and transposition operations on the arrangement of the datavolume in memory during the various steps of time migration. Note thatthe sequence depicted in Table II also allows data to be optionallywritten to disk if it is desired to stop and later restart the analysis.

The data manipulation in the invention is designed to take advantage ofthe structure of the MPP. Initially, the data are in a (t, x, y) volume,Step 1. In contrast to the generally horizontal frequency slicesdepicted in FIG. 2 for depth migration, these slices are generallyvertical, having a constant value along the y-axis. The data volume issliced into y-planes, so that a group of processors has a t-x plane. Thedata are transposed in MPP memory, Step 2, and Fourier transformed togive a (k, t, y) volume, Step 3.

Next, a sequence of three transpositions, Steps 4, 5, and 6, allow aFourier transform from the y spatial domain to the k_(y) wavenumberdomain, Step 7. One additional transposition, Step 8, facilitates aFourier transform from temporal into frequency space (ω, k_(y), k_(x)),Step 9. The data are now in a position to be time-migrated, Step 10. Aswas noted above, the order in which the transpositions and the Fouriertransforms are performed takes advantage of the structure of the MPP andthe memory allocation in it.

After the Step 10 time migration process, the transpositions andtransforms facilitate the return to the spatial domain [(τ, x, y)space]. These steps reverse the sequence of operations performed inSteps 2-8.

As indicated above in association with the depth migration embodiment ofthe present invention, the above time migration data manipulation schemeallows each PE to compute a direct solution of the wave equation in themigration step and to avoid the prior art's reliance on indirect,look-up table solutions in which, for example, separate operators wererequired for each frequency. The time migration embodiments of thepresent invention thereby attain more accurate, more cost-effectiveresults than have previously been possible.

                  TABLE II                                                        ______________________________________                                        Parallelization of Time Migration Data Manipulation Tasks                     Processors are in the 1-2 Plane.                                                                   Axis                                                     Step       Task      1          2   3                                         ______________________________________                                        1          Start     t          x   y                                         2          Transpose x          t   y                                         3          FFT       k.sub.x    t   y                                         4          Transpose*                                                                              t          k.sub.x                                                                           y                                         5          **Transpose                                                                             t          y   k.sub.x                                   6          Transpose y          t   k.sub.x                                   7          FFT       k.sub.y    t   k.sub.x                                   8          Transpose t          k.sub.y                                                                           k.sub.x                                   9          FFT       ω    k.sub.y                                                                           k.sub.x                                   10         Migrate   τ      k.sub.y                                                                           k.sub.x                                   11         Transpose k.sub.y    τ                                                                             k.sub.x                                   12         Inverse FFT                                                                             y          τ                                                                             k.sub.x                                   13         Transpose*                                                                              τ      y   k.sub.x                                   14         **Transpose                                                                             τ      k.sub.x                                                                           y                                         15         Transpose k.sub.x    ω                                                                           y                                         16         Inverse FFT                                                                             X          τ                                                                             y                                         17         Transpose τ      x   y                                         ______________________________________                                         *Optional write to disk followed by                                           **Optional read from disk                                                

Processing Overview

Data Flow--Generally

The overall flow of data is given in FIGS. 7A-7E. These Figures broadlyfocus on transformation of the x coordinate to the k_(x) domain (FIG.7A, steps 156-164), transformation of they coordinate to the k_(y),domain (FIG. 7B, steps 172-182), and transformation of the t coordinateto the frequency domain prior to imaging, along with an optional finalpre-imaging filter (FIG. 7B, step 182 and FIG. 7D, steps 220-224). Theimaging sequences are outlined on the remainder of FIG. 7D and on FIG.7E, and the post-imaging procedures required to transform the data backto the migrated time, generally referred to in the art as τ, and x and ydomains are shown on FIG. 7B (steps 184-190) and FIG. 7C.

The data are initially in a (t, x, y) cube, 150. The first step is todivide the data into y-slices, e.g., into planes in the t-x domain, 152.In contrast to the generally horizontal frequency slices depicted inFIG. 2 for depth migration, these slices are generally vertical, havinga constant value along the y-axis. An initial allocation of these slicesto processor groups is performed, 154. As above, it is to be understoodthat groups can consist of one or more processors. Next, the data arechecked to see if any slices remain unprocessed, 156. Initially, ofcourse, that would be the case and the next slice of t-x data is readinto the next group of processors, 158.

An FFT is applied in the x direction to give y-slices of t-k_(x) data,160. As discussed above, this FFT is most efficiently implemented on anMPP if each processor performing the calculation contains all data alongthe x direction for a specific y-slice. Therefore, the FFT requires atransposition of the t-x data into x-t data, application of the FFT inthe x-direction and a transposition back to give t-k_(x) data. Dependingon whether the PE can retain the data within its memory and continueprocessing, 162, processing either returns to 156, or writes the dataout to disk, 164, prior to going back to 156.

After all the y-slices of data have been processed into the k_(x)domain, transformation of the y coordinate to the k_(y) domain isperformed in an analogous manner (FIG. 7B, steps 172-182). However, thisportion of the processing requires that the (t, k_(x), y) data besubdivided into generally vertical k_(x) -slices, 172. As with theinitial y-slicing, the number of k_(x) -slices, and the allocation ofthose slices to processor groups, 172, will depend on the size of thedataset (i.e. the number of samples in each of the three directions t,k_(x) and y), the number of available processors, and the amount ofmemory available per processor, as will be understood to those skilledin the art.

If the data are in memory, 174, then they are transposed into k_(x)-slices, 176. Step 178 checks to see if all k_(x) -slices have beenprocessed, and if so, the final processing sequence in FIG. 7C isinitiated. If any of the k_(x) -slices is unprocessed, then the nextgroup of k_(x) -slices are read into memory from disk if necessary, 180.During this read from disk, the transpose to k_(x) slices whichpreviously occurred in step 176 for data in memory occurs. At this stageof the process, the data are in (t, y, k_(x)) form. Step 182 involvesthe final transpositions and transforms to (t, k_(y), k_(x)), and thenby FFT to (ω, k_(y), k_(x)) space.

The data are then processed following the scheme in FIGS. 7D and 7E, toproduce a migrated image in (τ, k_(y), k_(x)) space, where τ is themigrated time. This migration process is given below in the discussionof FIG. 7D and 7E. The procedure repeats via the additional sliceprocessing sequence of steps 188 and 178, with, if necessary, datawritten to disk transposing y and k_(x), 190. Note that steps 186 and190 perform the transform and transpositions in steps 11, 12, 13, and 14of Table II.

FIG. 7C outlines the steps involved in transformation from (ω, y, k_(x))space to (τ, x, y) space, and is similar to and the reverse of the stepsinvolved in the transpositions and transformations leading up to themigration step. Once again, the data volume is partitioned intoy-slices, 200, and the data are divided among processor groups, 202. Ifthe data are in memory, 204, they are transposed to give y-slices, 206.The slices are checked to see if any more y-slices remain, 208; if notthe process is terminated. If, however, there are additional slices tobe processed, they are read in to the processor groups from disk asy-slices if necessary, 210. The data are then transformed into the (τ,x, y) space, 212. This requires a transposition, an inverse FFT, andanother transposition to implement efficiently on an MPP. The migratedtime slices are written out to disk, 214, as the end product of themigration process.

Data Flow--Migration

FIGS. 7D and 7E show the migration process. The (ω, k_(y), k_(x)) datavolume from 182, FIG. 7B, is optionally filtered in the ω-domain, andthe velocity model divided into layers of constant velocity, 222. Notethat these constant velocity layers are less complicated than thevelocity model of the depth migration embodiment discussed above, andtherefore the communication of the velocities to the processors will notgenerally require the broadcast technique employed in depth migration,although that technique may be employed. The desired migrated timedomain τ is also divided into layers, 224. If the phase shift algorithmin the second embodiment is not to be used, 226, processing steps areshown in FIG. 7E. If the phase shift algorithm is to be used, the dataare checked to see if there are any remaining τ-layers to be processed,228. If not, the frequency domain migration is complete and the data aresent back to 184 in FIG. 7B. If additional τ-layers remain to beprocessed, the velocity for this layer is obtained; if it is differentfrom the velocity for the previous τ-layer, the frequency dependentphase shifts are computed, 230 and applied, 232, to downward continuefrom migrated time τ to time τ+Δτ. The data are summed over so to forman image at τ+Δτ, 234. This is repeated from 228 for additional Δτintervals.

The third embodiment performs a recursive Stolt migration, shown in FIG.7E. A full discussion of the conventional Stolt algorithm is given inStolt, R. H., Migration by Fourier Transforms, 43 Geophysics 23 (1978).Other methods are known in prior art and the disclosure herein is notmeant to be restricted to the recursive Stolt algorithm. The thirdembodiment here assumes that the velocity variation in the verticaldirection may be approximated by layers within which the velocity isconstant. The recursive Stolt method is also discussed in the StatutoryInvention Registration H 482 of Berryhill, Gonzalez and Kim issued in1988. However, the third embodiment is different from the inventiondisclosed by Berryhill et. al. in the method used for time migration.While the effect of FFT steps in the present invention is similar tothat taught by Berryhill et. al., their invention does not have thematrix transpositions between successive Fourier transforms that isinvolved in the present method.

A check is made for any remaining velocity layers, 236. If no velocitylayers remain to be processed, the frequency domain migration iscomplete and processing goes back to 184. If additional velocity layersremain to be processed, the velocity for this V-layer is obtained, 240.A Stolt migration is performed for all τs in this velocity layer andbelow, 242. There is no explicit downward continuation at this point.The τs in this velocity layer are saved, 244, discarding deeper migratedtimes. Next, the unmigrated data from the top of the layer are downwardcontinued to the bottom of the layer using the layer velocity. Thisrequires computation of a phase shift, 246, followed by a downwardcontinuation, 248. This is repeated for additional τ-layers, going backto 184 in FIG. 7B at the end of the completion of thefrequency-wavenumber migration.

Imaging Calculations

The phase shift method of Gazdag, J., Wave-Equation Migration with thePhase Shift Method, 43 Geophysics 1342, (1978) is the basis for the timemigration process in the second embodiment. The zero offset seismicsection, p (x, t, τ=0) for the 2-D case, may be considered as awavefield measured at the surface of the earth. The variables x, t, andτ are the horizontal position, the two-way travel time, and the two-waymigrated traveltime, respectively. The migration is performed by usingthe expression for downward continuation in the frequency domain

    P(k.sub.x, ω, τ+Δτ)=P(k.sub.x, ω, τ) exp (-i k.sub.z v Δτ/2)                                 (2),

where

    k.sub.z =(4 ω.sup.2 /v.sup.2 -k.sub.x.sup.2).sup.1/2 (3),

with k_(x) being the horizontal wavenumber and v the velocity in theinterval Δτ.

This invention uses equation (2) in a layer by layer extrapolation,replacing k_(x) ² in equation (3) by (k_(x) ² +k_(y) ²) for the 3-Dequivalent of Gazdag's derivation. The time-migrated image is obtainedas p(x, t=0, τ), which is referred to by those skilled in the art as theimaging condition. In the frequency domain, this is obtained bysummation of the various frequency components in the solution ofequation (2), at 234 of FIG. 7D. As would be known to those familiarwith the art, alternate solutions to the wave equation could be used,rather than equation (2). The present implementation of equation (2) isnot meant to be a limitation on other means for wavefield extrapolation.

The third embodiment uses a recursive frequency-wavenumber migrationprocess based on the so called Stolt migration. Stolt migration uses theresult that the downward continued P(k_(x), τ, ω) in terms of theunmigrated wavefield P(k_(x), 0, ω) and an extrapolation operator E(τ)is given by

    P(k.sub.x, τ, ω)=E(τ)P(k.sub.x, 0, ω)  (4).

The extrapolation operator is given by

    E(τ)=exp[-i(ω.sup.2 -W.sup.2).sup.1/2 τ]h(ω.sup.2 -W.sup.2)                                                 (5),

where h(. . . ) is the unit step function and

    W=k.sub.x v/2                                              (6),

v being the velocity in the layer.

The time migrated image is obtained as the inverse Fourier transform ofΩ/(Ω² +W²)^(1/2) P(k_(x), τ, Ω) where

    Ω=ω(1-W.sup.2 /ω.sup.2).sup.1/2          (7).

The recursive scheme first does the migration with a velocityappropriate for the first layer, discarding migrated times below thislayer. This is followed by a downward continuation of the unmigratedwavefield to the bottom of the first layer. The downward continuation isdone by the phase shift method given in the second embodiment. Thisproduces data that would have been recorded at the bottom of the firstlayer. The Stolt migration is then performed with the velocity for thenext layer and the process of downward continuation and migrationrepeated until the desired migrated time has been reached.

The recursive migration method is significantly faster than priormethods of migration because the migration is performed by making use ofan FFT, as opposed to a convolution, and the migration and the downwardcontinuation are performed over relatively large intervals in which thevelocity is constant. The second embodiment requires a summation of theextrapolated wavefield over all frequencies to give a migrated image inthe frequency domain. The third embodiment gives the time-migrated imagefor all τ within the velocity layer by an inverse Fourier transform. Thephase shift used in the downward continuation is the same in the secondand third embodiments. When the velocity layer thickness is equal to theτ-layer thickness, the two methods are theoretically equivalent,although computational efficiency will vary depending on implementationconsiderations, such as the FFT calculation routine being employed.

Depth Migration Example

The performance of the present invention, as compared to the prior art,is demonstrated in FIG. 8, which illustrates the reduction in clock timefor a 3-D post stack depth migration using the present methodimplemented on a T3D MPP compared with that for a prior art methodimplemented on a Y-MP. FIG. 8 shows the clock time for a relativelysmall 3-D survey for which a post-stack depth migration was implementedon a 4 processor Cray Y-MP. Even for the relatively small data volume of1×10¹¹, more than 60 hours were required to perform the migration. Incomparison, a data volume almost five times as large can be processed ona 128 processor T3D in less than 20 hours. A data volume of 4×10¹¹, forexample, corresponds to 100 frequencies for a survey size of 20 km.×20km. with a sampling every 10 meters, and 1000 depths every 5 meters. Theresult is dramatically lower seismic data processing costs.

It should be understood that the invention is not to be unduly limitedto the foregoing that has been set forth for illustrative purposes.Various modifications and alternatives will be apparent to those skilledin the art without departing from the true scope of the invention.

                  APPENDIX                                                        ______________________________________                                        C TITLE - Finite difference coefficient generator (TCDGEN.sub.-- Y)           C FUNCTION - This subroutine generates the linear systems that result         from the                                                                      C   2D 65-degree finite difference approximation to each of the                   crosslines.                                                               C   This routine handles both the case where a frequency plane spans              several                                                                   C   PEs and the case where it is restricted to one PE.                        C LANGUAGE - FORTRAN                                                          C LINKAGE - FORTRAN CALL STATEMENT                                            C   CALL TCDGEN.sub.-- Y (IER, NX, NY, OMEGA, VELOCITYD2,                         CP,                                                                       C    DZ, DY, OFFD, DIAG, RHS, INDEX, MYPE, NP, NLIST,                              MSGID)                                                                   C   PARM    DESCRIPTION I/O                                                   C   IER Error return   I I/O                                                  C     = 0 No error                                                            C     = 1 MYPE is not an element of NLIST                                     C     = 2 NY < 1                                                              C     = 3 NP < 1                                                              C     = 4 Error explained with a write statement                              C   NX Number of bins each line.   F I                                        C   NY Number of gridlines on this PE   F I                                   C   OMEGA Circular frequency for this slice.  E I                             C   VELOCITYD2 is the velocity/2 field, size (NX,NY). E*4 I                   C   CP Complex freq. slice, size (NX,NY)  C*8 I                               C   DZ Vertical step size.    E I                                             C   DY Spacing between crossline traces.   E I                                C   OFFD Complex array for off-diagonal of matrices C*8 O                     C     of size NY                                                              C   DIAG Complex array for diagonal of matrices C*8 O                         C     of size NY.                                                             C   RHS Complex array for right hand side of  C*8 O                           C     linear systems. Dimensioned (NY),                                       C   INDEX The gridline to be processed   I I                                  C   MYPE My processor id   I I                                                C   NP  The number of PEs in this group   I I                                 C   NLIST(NP)                                                                 C     A list of the PEs in this group   I I                                   C   MSGID The next valid message identifier   I I/O                           C EXT. REF.- PRTERR, PVMFPRECV, PVMFPSEND                                       SUBROUTINE TCDGEN.sub.-- Y (IER, NX, NY, OMEGA,                               VELOCITYD2, CP, DZ,                                                           &   DY, OFFD, DIAG, RHS, INDEX,                                               &   MYPE, NP, NLIST, MSGID)                                                   IMPLICIT NONE                                                                 INTEGER IER, NX, NY, INDEX, MYPE, NP, NLIST(NP), MSGID                        REAL OMEGA, DZ, DY                                                            REAL*4 VELOCITYD2(NX,NY)                                                      COMPLEX*8 CPJ, CPJP1, TEMP1, TEMP2                                            COMPLEX*8 CP(NX,NY), A1, A2, RHS(NY), OFFD(NY),                               DIAG(NY), C1                                                                  COMPLEX CP0, CPNYP1, CPFIRST, CPLAST                                          CHARACTER*500 ERR.sub.-- MSG                                                  INTEGER I,J, MY.sub.-- POS, INFO, IWHO, IWHAT, ITAG                           REAL*4 EPS, GAMMA, ALPHA, BETA, ONE, TWO, C2                                  REAL B1, B2                                                                   INCLUDE `fpvm3.h`                                                             IER = 0                                                                       EPS = 1.0E-14                                                                 GAMMA = 0.14                                                                  ALPHA = 0.4782406                                                             BETA = 0.376369527                                                            ONE = 1.                                                                      TWO = 2.                                                                    C Check for nonsensical NP.                                                     IF(NP.LT.1) THEN                                                               IER = 3                                                                       CALL PRTERR (IER, MYPE,                                                      &  `Routine TCDGEN.sub.-- Y has been accessed with NP less                    than 1.`,                                                                     &  `TCDGEN.sub.-- Y`)                                                          GO TO 990                                                                    END IF                                                                      C Check for nonsensical NY.                                                     IF(NY.LT.1) THEN                                                               IER = 2                                                                       CALL PRTERR (IER, MYPE,                                                      &  `Routine TCDGEN.sub.-- Y has been accessed with NY less                    than 1.`,                                                                     &  `TCDGEN.sub.-- Y`)                                                          GO TO 990                                                                    END IF                                                                      C Find out which PE I am on this list.                                          IER = 1                                                                       DO I = 1,NP                                                                    IF(MYPE.EQ.NLIST(I)) THEN                                                     MY.sub.-- POS = 1                                                             IER = 0                                                                       END IF                                                                       END DO                                                                        IF(IER.NE.0) THEN                                                              CALL PRTERR (IER, MYPE,                                                      &  `MY processor is not in the group. This is in routine                      &  TCDGEN.sub.-- Y.`, `TCDGEN.sub.-- Y`)                                       GO TO 990                                                                    END IF                                                                      C Generate coefficient matrices and rhs for tridiagonal systems.              C Along Y - direction.                                                        C OFFD is the upper and lower diagonal of the tridiagonal MATRIX              C on the left hand side.                                                      C DIAG is the diagonal of the coefficient matrix.                             C DIAG is the main diagonaL.                                                    C1 = CMPLX (0.0, 1.0) * DZ * ALPHA / (2.0 * DY * DY *                         OMEGA)                                                                        C2 = BETA / (DY * DY * OMEGA * OMEGA)                                         DO J = 1, NY                                                                   OFFD(J) = VELOCITYD2(INDEX,J) *                                               (C2*VELOCITYD2(INDEX,J)                                                      &  -C1) + GAMMA                                                                DIAG(J) = ONE - TWO * OFFD(J)                                                END DO                                                                      C At the seams, exchange info for RHS and, possibly, for the                  C boundary conditions:                                                        C Shovel values to the left:                                                    IF(MY.sub.-- POS.GT.1) THEN                                                    CPFIRST = CP(INDEX,1)                                                         CALL PVMFPSEND(NLIST(MY.sub.-- POS-1),MSGID,CPFIRST,1,                        COMPLEX16,INFO)                                                              IF (INFO .LT. 0) THEN                                                          WRITE(ERR.sub.-- MSG`("PVMFPSEND has returned an error",                     &  "code of",19,". In routine TCDGEN.sub.-- Y.")`) INFO                        CALL PRTERR (INFO, MYPE, ERR.sub.-- MSG `PVMFPSEND`)                          IER = 4                                                                       GOTO 990                                                                      END IF                                                                       END IF                                                                      C Shovel values to the right:                                                   IF(MY.sub.-- POS.LT.NP) THEN                                                   CPLAST = CP(INDEX,NY)                                                         CALL PVMFPSEND(NLIST(MY.sub.-- POS+1),MSGID+1,CPLAST,                         1,COMPLEX16,                                                                 &  INFO)                                                                       IF (INFO .LT. 0) THEN                                                         WRITE(ERR.sub.-- MSG,`("PVMFPSEND has returned an error",                    & "code of",19,". In routine TCDGEN.sub.-- Y.")`)INFO                          CALL PRTERR (INFO, MYPE, ERR.sub.-- MSG `PVMFPSEND`)                          IER = 4                                                                       GOTO 990                                                                      END IF                                                                       END IF                                                                      C Catch values from the neighbors, form the remainder of the RHS, and         C set up the boundary conditions:                                             C Catch values from the right:                                                   IF(MY.sub.-- POS.LT.NP) THEN                                                  CALL PVMFPRECV(-1,MSGID,CPNYP1,1,COMPLEX16,                                  & IWHO,IWHAT,ITAG,INFO)                                                        IF (INFO .LT. 0) THEN                                                          WRITE(ERR.sub.-- MSG,`("PVMFPRECV has returned an error",                 & "code of",19,". In routine TCDGEN.sub.-- Y.")`) INFO                           CALL PRTERR (INFO, MYPE, ERR.sub.-- MSG, `PVMFPRECV`)                         IER = 4                                                                       GOTO 990                                                                      END IF                                                                     C Handle the left side of the seam and set up left boundary condition as      needed:                                                                          IF(MY.sub.-- POS.NE.1) THEN                                                   RHS(NY) = CONJG (OFFD(NY)) * (CPNYP1 + CP(INDEX,                              NY-1))                                                                       &  + (1.0-2.0 * CONJG (OFFD(NY))) * CP(INDEX,NY)                               ELSE                                                                          IF(NY.EQ.1) THEN                                                               A1 = 2.0 * CP(INDEX,1) * CONJG (CPNYP1)                                       B1 = CP(INDEX,1) * CONJG (CP(INDEX,1))                                      &  + CPNYP1 * CONJG (CPNYP1)                                                    IF(ABS (B1) .LE. EPS) THEN                                                    A1 = CMPLX (1.0, 0.0)                                                        ELSE                                                                          A1 = A1/B1                                                                    END IF                                                                        IF (AIMAG (A1) .LT. 0.0)A1 = CONJG(A1)                                        DIAG(I) = OFFD(1) * A1 + DIAG(1)                                              RHS(1) = CONJG (OFFD(1)) * CPNYP1                                            &  + ((1.0-2.0 * CONJG (OFFD(1)))                                             &  + CONJG (OFFD(1)) * A1) * CP(INDEX,1)                                       ELSE                                                                          RHS(NY) = CONJG (OFFD(NY)) * (CPNYP1 + CP(INDEX,                              NY-1))                                                                       &  + (1.0-2.0 * CONJG (OFFD(NY))) * CP(INDEX,NY)                               A1 = 2.0 * CP(INDEX,1) * CONJG (CP(INDEX,2))                                  B1 = CP(INDEX,1) * CONJG (CP(INDEX,1))                                       &  + CP(INDEX,2) * CONJG (CP(INDEX,2))                                         IF(ABS (B1) .LE. EPS) THEN                                                     A1 = CMPLX(1.0, 0.0)                                                         ELSE                                                                           A1 = A1 / B1                                                                 END IF                                                                        IF (AIMAG (A1) .LT. 0.0) A1 = CONJG (A1)                                      DIAG(1) = OFFD(1) * A1 + DIAG(1)                                              RHS(1) = CONJG (OFFD(1)) * CP(INDEX,2)                                       &  + ((1.0-2.0 * CONJG (OFFD(1)))                                             &  + CONJG (OFFD(1)) * A1) * CP(INDEX,1)                                       END IF                                                                        END IF                                                                       END IF                                                                      C Catch values from the left:                                                   IF(MY.sub.-- POS.GT.1) THEN                                                    CALL PVMFPRECV(-1,MSGID+1,CP0,1,COMPLEX16,IWHO,                               IWHAT,ITAG,INFO)                                                              IF (INFO .LT. 0) THEN                                                         WRITE(ERR.sub.-- MSG,`("PVMFPRECV has returned an error",                    & "code of",19,".In routine TCDGEN.sub.-- Y.")`) INFO                          CALL PRTERR (INFO, MYPE, ERR.sub.-- MSG, `PVMFPRECV`)                         IER = 4                                                                       GOTO 990                                                                      END IF                                                                     C Handle the right side of the seam and set up right boundary condition       as needed:                                                                       IF(MY.sub.-- POS.NE.NP) THEN                                                   RHS(1) = CONJG (OFFD(1)) * (CP(INDEX,2) + CP0)                              &  + (1.0-2.0 * CONJG (OFFD(1))) * CP(INDEX,1)                                 ELSE                                                                           IF(NY.EQ.1) THEN                                                              A2 = 2.0 * CP(INDEX,NY) * CONJG (CP0)                                         B2 = CP(INDEX,NY) * CONJG (CP(INDEX,NY))                                    &  + CP0 * CONJG (CP0)                                                         IF(ABS (B2) .LE. EPS) THEN                                                     A2 = CMPLX (1.0, 0.0)                                                        ELSE                                                                           A2 = A2/B2                                                                   END IF                                                                        IF(AIMAG (A2) .LT. 0.0) A2 = CONJG (A2)                                       DIAG(NY) = OFFD(NY) * A2 + DIAG(NY)                                           RHS(NY) = CONJG (OFFD(NY)) * CP0                                             &  + ((1.0-2.0 * CONJG (OFFD(NY)))                                            &  + CONJG (OFFD(NY)) * A2) * CP(INDEX,NY)                                     ELSE                                                                          RHS(1) = CONJG (OFFD(1)) * (CP(INDEX,2) + CP0)                               &  + (1.0-2.0 * CONJG (OFFD(1))) * CP(INDEX,1)                                 A2 = 2.0 * CP(INDEX,NY) * CONJG (CP(INDEX,NY-1))                              B2 = CP(INDEX,NY) * CONJG (CP(INDEX,NY))                                     &  + CP(INDEX,NY-1) * CONJG (CP(INDEX,NY-1))                                   IF(ABS (B2) .LE. EPS) THEN                                                     A2 = CMPLX (1.0, 0.0)                                                        ELSE                                                                           A2 = A2 / B2                                                                 END IF                                                                         IF(AIMAG (A2) .LT. 0.0) A2 = CONJG (A2)                                       DIAG(NY) = OFFD(NY) * A2 + DIAG(NY)                                           RHS(NY) = CONJG (OFFD(NY)) * CP(INDEX,NY-1)                                 &  + ((1.0-2.0 * CONJG (OFFD(NY)))                                            &  + CONJG (OFFD(NY)) * A2) * CP(INDEX,NY)                                      END IF                                                                       END IF                                                                       END IF                                                                      C Finish the RHS: (The compiler needs to ignore this loop when                NY < 2.)                                                                        CPJ = CP(INDEX,2)                                                             TEMP1 = CPJ - CP(INDEX,1)                                                     DO J = 2, NY - 1                                                               CPJP1 = CP(INDEX,J+1)                                                         TEMP2 = CPJP1 - CPJ                                                           RSH(J) = CPJ + CONJG(OFFD(J))*(TEMP2 - TEMP1)                                 CPJ = CPJP1                                                                   TEMP1 = TEMP2                                                                END DO                                                                      990 MSGID = MSGID + 2                                                           RETURN                                                                        END                                                                         ______________________________________                                    

We claim:
 1. A frequency domain method of processing seismic data, saidseismic data corresponding to a subsurface region of the earth, saidprocessing performed on a computer system having multiple processingelements, said method comprising the steps of:a) specifying slices insaid data; b) assigning each of at least a plurality of said processingelements at least one partition on at least one of said slices, whereinsaid partition assignment for each of said plurality of said processingelements is independent of said partition assignment of each other ofsaid plurality of said processing elements; c) precomputing a velocitymodel corresponding to said subsurface region; d) in each of saidplurality of said processing elements, performing a direct migration ofsaid partition in said element using said velocity model, said migrationperformed in the frequency domain, each of said plurality of saidprocessing elements performing said migration independent of each otherof said processing elements, said migration facilitated by Fouriertransformation of each of two spatial parameters characterizing saidseismic data into the wavenumber domain, each of said transformationsperformed by each of said plurality of processing elements independentof each other of said plurality of processing elements; e) generating amap of reflectors in said seismic data using said migrated partitions.2. The method of claim 1, whereina) said seismic data is Fouriertransformed from the time domain to the frequency domain and said slicesare frequency slices in said transformed data; b) said velocity modelextends to a maximum depth of interest in said subsurface region and issubdivided into a family of velocity slices; c) for each of saidplurality of processing elements, broadcasting a first of said assignedpartitions to said element; d) said migration involves(i) broadcasting afirst of said family of pre-computed velocity slices to each of saidplurality of processing elements; (ii) downward continuing each saidpartition in each said element using said broadcast velocity slice,wherein said downward continuation involves a split wave equationsolution technique in which a first of said two spatial parameters issolved for each said partition by said assigned processing elementindependent both of each other said partition and of each other saidprocessing element; and (iii) repeating said velocity slice broadcastingand said downward continuing procedures until said maximum depth ofinterest has been attained; e) repeating said partition broadcasting andsaid migration until all said partitions on all said slices have beenmigrated; and f) said map is a reflector depth map.
 3. The method ofclaim 2 wherein, for a second of said two spatial parameters, saiddownward continuation involves for each said frequency slice a burn atboth ends dual solution technique which allows simultaneous calculationby at least two processing elements for the solution of said split waveequation for said second spatial parameter for said slice.
 4. The methodof claim 2, wherein said downward continuation also involves an errorcorrection calculation in the wavenumber domain.
 5. The method of claim4, wherein for each frequency slice said independent transforms arefacilitated by at least one transpose of data for one of said spatialparameters between said processing elements which are assigned saidpartitions on said slice.
 6. The method of claim 2, wherein saiddownward continuation also involves a wraparound filter calculation inthe wavenumber domain.
 7. The method of claim 6, wherein for eachfrequency slice said independent transforms are facilitated by at leastone transpose of data for one of said spatial parameters between saidprocessing elements which are assigned said partitions on said slice. 8.The method of claim 1, whereina) said specified slices are time slicesin said seismic data; b) said velocity model involves a family ofconstant velocity layers, said model extending to a maximum depth ofinterest in said subsurface region, said family of layers communicatedto each of said plurality of said processing elements; c) for each ofsaid plurality of processing elements, broadcasting a first of saidassigned partitions to said element; d) said migration involves, in eachof said plurality of said processing elements,i) Fourier transformingsaid broadcast partitions from the time domain to the frequency domain,and ii) time migrating said transformed partition using said family ofvelocity layers until said maximum depth of interest has been attained;e) repeating said migration until all said partitions on all said sliceshave been migrated; and f) said map is a reflector time map.
 9. Themethod of claim 8, wherein said independent transforms are facilitatedby at least one transpose of data between said processing elements whichare assigned said partitions on said slice.
 10. The method of claim 8,wherein said direct migration employs a phase shift downwardcontinuation.
 11. The method of claim 8, wherein said direct migrationemploys a recursive frequency--wavenumber downward continuation.